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I.  INTRODUCTION  AND  SUMMARY 


1 . 1 BACKGROUND 

The  objective  of  the  Systems,  Science  and  Software  (S3) 
research  program  is  to  examine  the  parameters  that  affect  the 
seismic  signals  from  underground  explosions  and  earthquakes. 
Attention  is  primarily  directed  to  those  features  of  the 
seismic  waveforms  that  discriminate  between  the  two  classes 
of  events  and  that  reliably  indicate  the  explosion  yield. 
Current  research  includes  empirical  studies  of  the  available 
data,  time  signal  analysis,  and  the  development  and  application 
of  theoretical  and  numerical  methods  for  modeling  earthquakes 
and  explosions.  Emphasis  is  on  the  last  of  these.  In  partic- 
ular, we  are  applying  techniques  for  theoretically  simulating 
the  far-field  signatures  of  simple  and  complex  seismic  sources. 

This  report  summarizes  the  work  done  during  the  sixth 
three-month  period  of  the  current  contract. 


1.2  SUMMARY  OF  RESEARCH  DURING  THIS  QUARTER 

Our  work  during  this  quarter  has  included  research  in 
a number  of  areas.  Research  projects  currently  underway  or 
completed  during  the  quarter  are  briefly  summarized  below. 

Source  Studies 

A.  Analysis  of  Small-Scale  Explosions  in  Grout 

A final  report,  SSS-R-77-3349 , "Seismic  Ground 
Motion  from  Free-Field  and  Underburied  Explosive  Sources," 
by  J.  T.  Cherry,  T.  G.  Barker,  S.  M.  Day  and  P.  L.  Coleman 
was  submitted  during  this  quarter.  The  abstract  is  repro- 
duced below: 

Small-scale  laboratory  experiments  were  conducted  and 
analyzed  to  study  the  effect  of  the  proximity  of  the  free 
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surface  on  the  seismic  ground  motions.  Two  classes  of  experi- 
ments were  done.  In  one  the  charges  were  far  from  the  free- 
surface  and  the  free-field  displacement-time  histories  were 
measured.  In  the  second  class  the  charges  were  near  the  sur- 
face and  were  either  fully  contained  or  formed  a crater.  The 
charges  were  0.25  g of  PETN  in  concrete  cylinders,  120  cm  in 
diameter  and  33  to  60  cm  thick.  In  all  experiments  displace- 
ments were  measured  30.5  cm  directly  below  the  charges.  The 
experiments  produced  consistent  and  repeatable  data.  A 
striking  feature  of  the  ground  motions  for  the  near  surface 
experiments  is  a large  long-period  negative  pulse  which  is 
present  whether  or  not  cratering  occurred. 

The  results  were  studied  by  comparing  to  numerical 
simulations  of  the  experiments  using  a Lagrangian  finite 
difference  program  and  published  properties  of  the  concrete 
and  PETN.  The  calculations  are  in  good  agreement  with  the 
laboratory  data,  providing  verification  of  both  the  constitu- 
tive models  and  the  methods. 

The  large,  long-period  negative  pulse  can  be  explained 
in  the  context  of  linear  elasticity.  It  is  due  to  the  near- 
field interaction  of  the  spherical  wave  front  with  the  free 
surface.  Therefore,  it  does  not  propagate  to  the  far-fiel  I 
and  is  not  important  for  teleseismic  magnitudes. 

B.  Analysis  of  Decoupling  Calculations  Done  by  Applied 
Theory,  Inc.  (ATI) 

A topical  report  SSS-R-78-3627,  "Analysis  of  Two 
Decoupled  Explosion  Simulations"  by  T.  C.  Bache  and  J.  F.  Masso 
was  submitted  in  draft  form  in  April.  The  abstract  is  repro- 
duced below: 

Two  axisymmetric  ground  motion  calculations  were  carried 
out  by  Applied  Theory,  Inc.  to  simulate  25  kt  decoupled  nuclear 
explosions  in  mined  cavities  in  salt.  One  cavity  was  spherical 
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with  a radius  of  66  meters  while  the  other  was  a 3/1  aspect 
ratio  ellipsoid  of  revolution.  Both  had  the  same  volume. 
Results  of  the  two  calculations  were  analyzed  to  determine 
the  character  of  the  teleseismic  body  and  surface  waves. 

The  spherically  symmetric  portion  of  the  field  is 
slightly  (20  to  25  percent)  smaller  for  the  spherical  cavity. 
Comparing  to  results  of  tamped  explosions,  the  decoupling 
factor  for  this  case  is  about  140  at  one  Hertz.  The  radiation 
from  the  ellipsoidal  cavity  is  substantially  perturbed  from 
spherical  symmetry;  the  maximum  S wave  amplitudes  are  nearly 
three  times  as  large  as  maximum  P wave  amplitudes  at  1.0  Hz. 
However,  theoretical  body  and  surface  wave  seismograms  indi- 
cate that  the  m^  and  Mg  values  are  not  substantially  different 
for  the  two  cavities. 


C.  Transmitting  Boundary  for  Finite  Difference 
Calculations 


We  have  recently  implemented  a transmitting  boundary 
condition  for  finite  different  stress  wave  calculations,  based 
on  paraxial  approximations  to  the  wave  equation  given  in,  for 
example,  Engquist  and  Majda  (1977)  and  Clayton  and  Engquist  (1977). 
Our  approach  has  been  to  first  separate  the  field  near  the 
boundary  into  its  irrotational  and  solenoidal  components , and 
then  to  apply  a second  order  paraxial  approximation  to  each  of 
these  components.  The  method  has  been  tested  for  two-dimensional 
plane  boundaries  and  has  been  found  to  be  quite  satisfactory 
for  transmitting  the  P wave  and  the  initial  S wave  arrival. 

However,  boundary  displacement  suffers  from  a spurious  drift 
at  late  time.  The  size  of  the  late-time  drift  is  sensitive 
to  details  of  the  difference  scheme  used.  After  much  experi- 
mentation, we  believe  we  have  minimized  its  contribution  for 
the  second-order  approximation.  Higher  order  approximations 
may  be  required  to  remove  its  effect. 
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The  method  has  not  yet  been  extended  to  curved  boundaries 
like  those  parallel  to  the  axis  in  axisymmetric  problems.  We 
plan  to  attempt  this  extension.  We  also  plan  to  test  the  method 
for  three-dimensional  grids. 

D.  Earthquake  Modeling  on  the  ILLIAC  IV  Computer 

A three-dimensional  finite  different  program  for 
modeling  earthquake  faulting  has  been  made  operational  on  the 
ILLIAC  IV  computer.  The  current  version  is  capable  of  handling 
a bilateral  fault  in  a homogeneous  medium.  Our  initial  set  of 
calculations  is  designed  to  study  the  importance  of  plastic 
material  behavior  on  the  fault  zone  and  determine  the  scaling 
of  the  radiation  field  with  the  fault  parameters. 

The  first  calculation  was  completed  during  this  quarter. 
This  is  for  a 3 km  x 3 km  bilateral  fault  in  a linear  elastic 
space.  The  output  tapes  were  sent  to  S3  where  they  were  read 
and  selected  quantities  were  plotted.  Considerable  effort  was 
required  to  develop  the  software  to  convert  and  plot  these  data 
tapes . 

Unfortunately,  detailed  analysis  of  the  output  revealed 
some  strange  behavior  in  a region  of  the  fault  plane  that  was 
eventually  traced  to  a coding  error.  The  calculation  is  being 
repeated.  Processing  of  output  data  from  this  and  subsequent 
calculations  will  be  much  easier  due  to  the  availability  of  the 
recently  developed  software. 

E.  Representation  Theorem  for  Analytic  Continuation 

of  Finite  Difference  Source  Calculations 

This  work  is  summarized  in  Section  1.4  and  described 
in  detail  in  Section  III. 

F.  One-Dimensional  Source  Calculations 

One-dimensional  source  calculations  done  at  S3  in 
the  past  were  initiated  with  a rock-gas  filled  cavity  at  the 
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vaporization  radius  (approximately  70  metric  tons  of  rock  are 
vaporized  per  kiloton  yield) . An  ideal  gas  equation  of  state 
was  used.  We  have  recently  developed  complete  state  equations 
for  tuff  and  granite,  valid  up  to  many  megabars,  so  the  cal- 
culation can  be  initiated  at  the  wall  of  the  original  cavity. 
These  equations  of  state  give  a more  realistic  cavity  driving 
pressure  over  the  entire  time  of  the  calculation.  Results 
from  the  old  and  new  methods  for  handling  the  cavity  region  are 
being  compared. 

With  a complete  equation  of  state  we  can  compute  the 
decoupling  effect  of  varying  the  initial  radius  of  an  air  filled 
cavity.  We  are  currently  beginning  a series  of  decoupling  cal- 
culations in  granite,  tuff  and  salt.  These  calculations  will 
also  use  a better  one-dimensional  treatment  of  the  overburden 
pressure  than  was  used  in  previous  calculations. 


Data  Analysis 

A.  Discrimination  Experiment 

This  work  is  summarized  in  Section  1.3  and  described 
in  detail  in  Section  II. 

B.  Development  of  Spectral  m^  and  Mg 

The  variable  frequency  magnitude  (VFK)  discriminant 
being  applied  in  the  discrimination  experiment  (see  Section  II) 
is  based  on  the  use  of  narrow  band  filter  magnitudes,  called 
n^(f).  These  magnitudes  reflect  the  amplitude  of  the  arriving 
energy  at  a particular  period  and  a particular  arrival  time. 

In  essence,  they  represent  the  spectral  energy  of  a specific 
phase  arrival.  We  have  long  recognized  that  measurements  of 
this  kind  might  provide  a more  convenient,  stable  and  reliable 
signal  amplitude  indicator  than  the  conventional  m^. 
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A set  of  HNME  recordings  of  eleven  Pahute  Mesa  recordings 
was  provided  in  digital  form  and  is  being  examined  to  test  the 
utility  of  the  rn^f).  A similar  value,  Mg(f),  is  defined  for 
the  surface  waves.  Single  frequency  values  are  no  improvement 
on  the  conventional  time  domain  magnitudes.  However,  m.  (f)  and 
Mg(f)  values  defined  by  averaging  the  m^f)  and  Mg(f)  over  a 
frequency  band  are  very  attractive  in  terms  of  such  qualities 
as  the  scatter  when  plotted  versus  yield.  We  are  continuing 
to  evaluate  this  measurement  and  will  report  the  results  to  VSC 
in  the  near  future. 


Surface  Wave  Studies 

A.  Analysis  of  Surface  Waves  from  the  French  Test  Site 
m the  Sahara 

This  work  is  summarized  in  Section  1.6  and  described 
in  detail  in  Section  V. 

B.  Surface  Wave  Amplitudes  of  NTS  Explosions  Recorded 
at  ALQ  and  TUC 

In  previously  published  work  (Bache,  Rodi  and 
Harkrider,  1978)  we  presented  crustal  models  for  the  NTS-ALQ 
and  NTS-TUC  paths  that  lead  t<5  synthetic  seismograms  that 
closely  match  the  observations.  The  amplitudes  of  the  syn- 
thetics can  then  be  used,  together  with  the  theory,  to  deduce 
the  long  period  levels  of  the  source  functions  for  NTS  explo- 
sions. This  work  is  nearly  complete  and  a report  describing 
the  results  is  in  preparation. 

C.  Theoretical  Study  of  the  Lq  Phase 

This  work  is  summarized  in  Section  1.5  and  described 
in  detail  in  Section  IV. 
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Body  Wave  Studies 

In  the  section  on  "Data  Analysis"  we  mentioned  the  set 
of  eleven  HNME  recordings  of  Pahute  Mesa  explosions  that  are 
being  used  to  test  the  development  of  newly  defined  magnitude 
measures.  These  data  are  also  being  analyzed  with  synthetic 
seismogram  methods  to  isolate  the  different  effects  contri- 
buting to  m^-log  yield  and  Mg-log  yield  relations.  An  important 
feature  of  the  body  wave  recordings  is  a depth-dependent  phase 
that  arrives  after  pP.  This  may  be  due  to  spall  slapdown, 
tectonic  release  or  some  other  cause . More  complex  synthetic 
seismograms  are  being  studied  to  determine  the  origin  and 
nature  of  this  phase.  We  need  to  learn  how  to  model  it  for 
body  waves  to  estimate  its  effect  on  the  surface  wave  recordings. 


1.3  SUMMARY  OF  SECTION  II:  "DISCRIMINATION  EXPERIMENT" 

During  this  reporting  period,  we  received  digital  short- 
and  long-period  body  and  surface  wave  data  for  nineteen  Eurasian 
events  recorded  at  a global  network  of  seismic  stations.  Our 
initial  efforts  at  event  discrimination  have  focused  on  the 
application  of  narrow-band  filtering  techniques  to  the  short- 
period  P waves  recorded  at  four  of  the  stations  in  the  network. 
These  four  stations  are  located  at  epicentral  distances  between 
16  and  98  degrees  from  the  events. 

Tentative  identifications  of  thirteen  of  the  larger  events 
are  as  follows:  earthquakes  - events  47  and  49;  explosion-like 
events  1,  14,  16,  17,  18,  19,  20,  21,  22,  33  and  53.  The  event 
locations  are  described  in  Figure  1 and  Table  1 of  Section  II. 
More  definitive  results  must  await  the  receipt  and  analysis  of 
additional  events  and  the  determination  of  the  extent  of  the 
earthquake  and  explosion  populations  in  the  m^ff)  plane  at  each 
of  the  reporting  stations. 
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1.4  SUMMARY  OF  SECTION  III:  "ANALYTIC  CONTINUATION  OF  THE 

ELASTIC  FIELD  FROM  A COMPLEX  SOURCE  IN  A HALFSPACE" 

An  important  problem  in  computing  theoretical  seismograms 
* for  complex  deterministic  source  models  is  the  development  of 

methods  for  linking  the  finite  difference  numerical  source  models 
with  efficient  analytical  techniques  for  propagating  elastic 
waves  in  layered  earth  models.  When  the  source  is  in  a homo- 
geneous whole  space,  accurate  techniques  are  available  and  have 
been  extensively  used  at  S3.  However,  the  problem  is  much 
harder  when  there  are  material  interfaces  in  the  source  region. 

In  particular,  there  is  a great  need  for  an  analytic  continua- 
tion method  for  source  calculations  including  the  free  surface 
in  the  source  region. 

In  Section  III  we  outline  mathematically  an  analytic 
continuation  method  for  a source  in  a layered  half space.  The 
method  is  based  on  an  elastodynamic  representation  theorem  and 
requires  integration  of  a series  of  terms  involving  products  of 
the  source  generated  displacements  and  tractions  with  Green's 
functions  representing  the  elastic  response  of  the  medium.  The 
method  is  now  being  programmed  for  computing  the  body  waves  at 
arbitrary  range  in  a plane-layered  earth  model.  The  extension 

to  surface  waves  is  straightforward. 

C 

1.5  SUMMARY  OF  SECTION  IV:  "THEORETICAL  COMPUTATION  OF  Lg" 

r A brief  study  of  the  excitation  of  the  regional  phase  Lg 

was  conducted  at  VSC  request  using  a continental  earth  model 
provided  by  VSC.  Theoretical  Lg  phases  are  associated  with 
stationary  phases  of  the  higher  mode  Rayleigh  and  Love  waves 
t which  occur  at  about  3.5  km/sec  in  continental  structures. 

Synthetic  seismograms  illustrating  the  Lg  phase  for  several 
source  depths  and  ranges  are  presented  in  Section  IV.  Unfor- 
tunately, the  structure  near  the  surface  is  much  too  simple  to 
t give  realistic  synthetic  seismograms . The  synthetics  are 


dominated  by  a nearly  undispersed  pulse  associated  with  the 


1.6  SUMMARY  OF  SECTION  V:  "ANALYSIS  OF  SURFACE  WAVES  FROM 

THE  FRENCH  TEST  SITE  IN  THE  SAHARA" 

Surface  wave  recordings  from  seismometer  stations  in  or 
near  Africa  were  used  to  study  the  French  underground  nuclear 
explosion  SAPHIRE  (27  February  1965)  with  the  objective  being 
to  determine  the  source  characteristics  of  this  event.  Our 
first  step  is  to  accurately  account  for  the  path  effects  by 
inverting  the  observed  dispersion  to  determine  average  crustal 
models  for  the  five  paths  studied.  Models  fitting  Rayleigh 
wave  phase  and  group  velocities  were  inverted  for  each  path. 

In  addition,  Love  wave  group  velocities  were  simultaneously 
inverted  for  two  of  the  five  paths.  The  fit  to  the  observed 
data  is  excellent  for  all  but  one  path  where  the  Love  and 
Rayleigh  wave  group  velocities  cannot  be  fit  simultaneously. 
The  inversion  models  are  smooth  and  give  profiles  for  the 
crust  and  upper  mantle  that  are  consistent  with  data  from 
other  sources.  These  models  are  preliminary  in  the  sense 
that  they  have  not  yet  been  integrated  to  give  a consistent 
picture  of  the  structure  in  northern  Africa.  This  is  the 
next  step  which  is  necessary  before  considering  our  models 
to  be  our  final  estimates  for  this  region. 

The  inversion  models  were  used  together  with  a simple 
source  to  compute  theoretical  seismograms  for  the  SAPHIRE 
event.  These  seismograms  agree  very  well  with  the  observa- 
tions at  several  of  the  stations,  but  there  is  significant 
disagreement  at  the  east  African  stations  in  Ethiopia  and 
Kenya.  The  problem  may  be  with  the  path  models,  the  source 
representation  or  the  anelastic  attenuation  (Q)  model  used. 
These  questions  require  further  study.  After  completing  our 
path  model  specification,  we  will  use  the  comparison  of  syn- 
thetic and  observed  seismograms  to  infer  the  long  period 
source  characteristics. 
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II.  DISCRIMINATION  EXPERIMENT 


2 . 1 INTRODUCTION 

Our  objective  in  the  discrimination  experiment  is  to 
analyze  seismic  waveforms  from  a large  population  of  Eurasian 
events  in  order  to  identify  these  events  as  either  earth- 
quakes or  explosions.  The  waveforms  that  we  are  concentra- 
ting on  are  short-period  P waves  recorded  by  a global  net- 
work of  seismograph  stations.  In  this  section  of  the  report 
we  will  summarize  the  work  that  has  been  performed  to  date  on 
this  experiment. 

2.2  DESCRIPTION  OF  DATA 

As  of  the  end  of  this  reporting  period  we  had  received 
multi-station  short-  and  long-period  digital  magnetic  tape 
data  for  nineteen  Eurasian  events.  The  event  locations  are 
shown  in  Figure  1 and  their  dates,  origin  times  and  epi- 
central  coordinates,  as  supplied  by  Teledyne  Geotech,  are 
listed  in  Table  1.  These  events  are  well  distributed  through- 
out the  Eurasian  study  region,  from  Kamchatka  and  the  Kuril 
Islands  in  the  east  to  the  Caspian  Sea  in  the  west.  The 
locations  of  eight  of  the  Eurasian  stations  contributing  data 
for  this  experiment  are  also  indicated  on  Figure  1. 

Upon  receipt,  the  magnetic  tapes  are  converted  to  a 
format  more  convenient  for  use  on  the  UNIVAC  1108  computer. 
Once  the  tape  conversion  process  has  been  completed,  plots 
of  the  short-period  P wave  seismograms  are  generated  for 
those  event-station  pairs  that  are  not  included  in  the  sets 
of  Xerox  copies  of  seismograms  supplied  by  Geotech.  These 
analog  recordings  are  necessary  for  culling  the  data  and 
defining  the  duration  of  signal  and  noise  windows.  Examples 
of  P-wave  seismograms  recorded  at  six  different  stations  for 
event  19  are  shown  in  Figure  2.  As  can  be  seen,  the  duration 
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Figure  1.  Map  of  Eurasia  showing  locations  of  events  and  several  of  the  stations 
included  in  the  discrimination  experiment. 


TABLE  1 


EVENT  LOCATION  INFORMATION 


Coordinates 


Event 

Number 

Yr 

Date 

Mo  Dy 

Origin  Time 

Hr :Mn:Sec 

Latitude 
( °N) 

Longitude 
( °E) 

1 

77 

07 

26 

16:59:59.9 

69.4 

90.4 

14 

77 

07 

30 

01:56:59.9 

49.7 

78.2 

16 

77 

08 

10 

22:00:00.7 

50.9 

111.0 

17 

77 

08 

17 

04:26:59.8 

49.8 

78.2 

18 

77 

08 

20 

22:00:00.6 

64.1 

99.8 

19 

77 

09 

01 

03:00:00 

73.0 

54.0 

20 

77 

09 

05 

03:03:00 

50.0 

78.9 

21 

77 

09 

10 

16:00:00 

57.0 

106.0 

22 

77 

09 

30 

06:59:00 

48.0 

48.0 

33 

77 

10 

09 

11:00:00 

73.0 

55.0 

38 

77 

10 

16 

15:02:49 

36.9 

71.5 

41 

77 

10 

13 

20:38:42 

38.1 

72.8 

47 

77 

10 

16 

21:05:35 

49.7 

155.1 

48 

77 

10 

19 

05:02:00 

36.3 

71.3 

49 

77 

10 

19 

21:20:37 

49.5 

155.4 

53 

77 

10 

29 

03:07:00 

49.0 

78.0 

55 

77 

10 

26 

05:38:52 

49.0 

155.8 

57 

77 

10 

26 

13:14:30 

51.5 

153.4 

59 

77 

10 

28 

21:15:11 

39.8 

71.9 

Figure  2.  Short-period  P-wave  seismograms  for  event  19  recorded  at  several  of 
the  stations  included  in  the  discrimination  experiment. 


of  noise  preceding  these  signals  is  quite  variable,  ranging 
from  about  20  seconds  at  TATO  (19050)  to  70  seconds  at  ILPA 
(19100) . In  addition,  these  analog  playouts  help  us  avoid 
all  or  portions  of  seismograms  that  contain  bad  data  (clipped 
signals,  digital  drop-outs,  etc.),  such  as  the  NORSAR  record- 
ing (19085)  at  the  bottom  of  Figure  2. 

2.3  ANALYSIS  PROCEDURE 

The  variable  frequency  magnitude  discriminant,  as 
proposed  by  Archambeau,  et  al. , (1974)  and  developed  by  Savino 
and  Archambeau  (1974) , is  the  technique  that  will  be  used 
to  identify  the  Eurasian  events  in  this  experiment.  The 
variable  frequency  magnitudes  are  outputs  of  the  MARS  com- 
puter program. 

The  central  feature  of  the  MARS  signal  analysis  pro- 
gram is  the  use  of  narrow  band  frequency  filters  to  break 
up  or  decompose  a time  series  consisting  of  signal  plus 
noise  into  a set  of  quasi-harmonic  modulated  "signals”. 

This  set  of  filtered  signals,  one  for  each  filter  center 
frequency,  can  then  be  used  to  determine  the  energy  arrival 
time  (the  group  time,  tg)  and  amplitude  of  the  signal  for 
each  center  frequency  by  analysis  of  the  time  modulation  of 
the  filter  outputs. 

Figure  3 summarizes  the  flow  of  operations  in  the 
MARS  program:  Figure  3a  provides  a verbal  outline;  Figure 
3b  presents  the  key  mathematical  operations  performed  in 
this  program. 

Seismic  data  are  read  into  the  program  in  the  form 
of  time  series  (Figures  3a  and  3b)  generally  of  about  500 
to  2000  points  in  length.  The  data  are  then  optionally 
detrended,  demeaned  and  tapered.  The  program  then  selects 
the  smallest  power  of  two  which  is  greater  than  the  number 
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OOP  OVER 
QUENCS  OF 
I SMOG RAMS 


READ  AND  EDI" 

? INPUT  DATA 

LIST : 

1 . CONTROL,  PLOT , S 
REPORT  OPTIONS 

2 . INPUT  DATA 

3.  N3F  PARAMETERS 


READ  DESIRED  EVENT, 
SINGLE  COMPONENT  TIME  SERIES 
? signal  and  NOISE 

R 


ADDITIONAL 

COMPONENTS 


SCALE,  DEMEAN,  DETREND 
AND  TAPER  TIME  SERIES 


PERFORM  FFT  OF  TIME  SERIES 


PLOT  SIGNAL 
SPECTRUM 


CORRECT  SIGNAL  SPECTRUM  FOR 
INSTRUMENT  RESPONSE 


PLOT  CORRECTED 
SIGNAL  SPECTRUM 


CALCULATE  NOISE  AND  SIGNAL 
POWER  SPECTRA:  SIGNAL-TO- 
NOISE  RATIO 

3?  OVER  DIFFERENT] 

ILTER  FREQUENCIES ^ (fv) 


PLOT  NOISE  AND 
SIGNAL  POWER; 
SIGNAL-TO- NOISE 
RATIO 


NARROW  BAND  FILTER  (NBF) 
CORRECTED  SIGNAL  SPECTRUM 


CALCULATE  QUADRATURE 
SPECTRUM  BY  HILBERT  TRANSFORM 


PERFORM  INVERSE  FFT 
CF  BOTH  FILTERED  AND 
QUADRATURE  SPECTRA 


USE  TIME- DOMAIN  FILTERED 
SIGNAL  AND  QUADRATURE  TO 
FORM  ENVELOPE 


CALCULATE  INSTANTANEOUS 
PHASE  AND  FREQUENCY 


PLOT  NBF  OUTPUT 
SIGNAL 


PLOT  NBF 
ENVELOPE 


SAVE  NBF  SIGNAL, 
ENVELOPE,  AND  PHASE 
FOR  POLARIZATION 
FILTERING 


PLOT  INSTANTANEOUS 
PHASE  AND  FREQUENCY 


LOCATE  PEAKS  OF  ENVELOPE 
AND  SORT  PEAKS  IN 
DESCENDING  ORDER 


SAVE  PEAK 
AMPLITUDE,  PHASE, 
AND  GROUP 
ARRIVAL  TIME 


COMPUTE  VARIABLE  FREQUENCY 
BODY  AND  SURFACE  WAVE 
MAGNITUDES 


SAVE  BODY  AND 
SURFACE  WAVE 
MAGNITUDES 


PRINT  SUMMARY  OF  NBF  RESULTS 


PLOT  GROUP  ARRIVAL 
TIMES  VERSUS 
FREQUENCY 


PLOT  SUM  OF  ALL  NBF 
ENVELOPES  VERSUS  TIME 


rigure  3a.  (continued) 


PERFORM  ANALYSIS  OF  COMBINED 
NARROW- BAND,  POLARIZATION,  AND 
DIRECTIONAL  FILTERING  TO 
IDENTIFY  SPECIFIC 
SEISMIC  PHASES 
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Figure  3a. 


(continued) 
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Figure  3b.  (continued) 
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Figure  3b.  (continued) 
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COMPUTE  POLARIZATION  FILTER  FACTORS  FOR  EACH  PAIR 

OF  COMPONENTS; 

E.G.,  THE  VERTICAL-RADIAL  (Z-R) 
POLARIZATION  FACTORS  ARE 


p£k*  (t)  ■ sin  A $2^  (t) 


CIRCULAR  FUNCTION 


P^k> (t)  = cos  A $ *k) (t) 
ER(k) (t)  = A^k) (t)/A^k) (t) 


LINEAR  FUNCTION 


ELLIPSE  RATIO 


ALF(k) (t)  = arctan  [A*k) (t)/A^k) (t) ] APPARENT  P-WAVE 

L EMERGENCE  ANGLE 

WHERE  A4^k)  = $<k)  ft)  - ^k)  (t) 


COMPUTE  POLARIZATION  AND  DIRECTIONAL  FILTER  DETECTORS 
FOR  SPECIFIC  TYPES  OF  WAVE  MOTION; 

E.G. , THE  P-WAVE  DETECTOR  IS 

P (k)  (t)  = cos'  1 A $^k)  (t) 

• cos  2 [arctan  (A^/aJ^  ) - c^]  , 
if  P^k)  (t)  > 0 

WHERE  Mt  AND  M-  ARE  INTEGERS  AND  ct..  IS  THE  APPARENT 

EMERGENCE  ANGLE.  1 


Figure  3b.  (continued) 


CONSTRUCT  POLARIZATION  FILTERED  OUTPUT  SIGNALS  FOR 
SPECIFIC  TYPES  OF  WAVES;  E.G.,  FOR  P-WAVES 

Z^k)  Ct)  = X^0  (t)  • P(i°  (t)  , 

WHEN  Xz  IS  THE  NBF  OUTPUT  SIGNAL  AND  P^  IS  THE 
P-WAVE  POLARIZATION  DETECTOR 


ft 


of  points  input  and  performs  a discrete  Fourier  transform 
using  the  Cooley-Tukey  algorithm.  Both  the  original  time 
series  and  the  spectrum  are  then  plotted. 

Referring  to  Figures  3a  and  3b,  the  signal  is  next 
filtered  in  the  frequency  domain  by  multiplication  by  a 
narrow  band  Gaussian-shaped  filter.  This  particular  filter 
form  is  selected  to  satisfy  two  goals:  (1)  minimum  width 
in  the  frequency  domain,  and  (2)  maximum  ripple  suppression 
in  the  time  domain.  As  a consequence  of  the  sampling  theorem, 
or  uncertainty  principle,  one  cannot  simultaneously  satisfy 
these  two  goals  to  arbitrary  precision.  The  filter  employed 
was  selected  for  its  optimal  time  and  frequency  domain  char- 
acteristics within  this  basic  limitation. 

Once  the  signal  has  been  narrow  band  filtered,  it  is 
corrected  for  the  appropriate  instrument  response;  the  filtered 
signal  transform  is  divided  by  the  instrument  transfer  function. 
The  resulting  complex  spectrum  is  then  inverse  Fourier  trans- 
formed into  the  time  domain,  to  produce  what  will  hereafter  be 
referred  to  as  the  filtered  signal. 

The  narrow  band  filtered  signal  will  appear  as  a quasi- 
sinusoidal  carrier  wave  contained  within  a smooth  envelope. 

The  next  step  in  the  program  (Figures  3a  and  3b)  is  to  con- 
struct the  envelope  function  by  means  of  the  Hilbert  trans- 
form. In  particular,  a quadrature  signal  is  formed  by  multi- 
plying the  transform  of  the  filtered  signal  by  -i  sgn  (w) 
and  then  inverting  this  to  the  time  domain  by  an  inverse  trans- 
formation. The  envelope  function  is  then  constructed  by  taking 
the  square  root  of  the  sum  of  the  squares  of  the  filtered 
signal  and  its  quadrature.  The  maximum  of  the  envelope  func- 
tion occurs  at  the  time  of  arrival  of  energy  at  the  center 
frequency  of  the  filter  and  the  amplitude  of  the  maximum  is 
proportional  to  the  spectral  amplitude  of  the  filtered  signal 
at  the  center  frequency  of  the  filter.  The  narrow  band 
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band  filtering  procedure  can  be  performed  on  a particular 
component  seismogram  at  a number  of  different  frequencies 
within  some  band  of  interest.  Correlation  of  the  resulting 
envelope  functions  indicates  the  arrival  times  of  the  various 
frequency  components. 

The  instantaneous  frequency  and  phase  are  computed 
from  the  filtered  signal,  and  its  quadrature  signal  as  well, 
and  are  stored  for  subsequent  use  in  polarization  filter- 
ing with  additional  components  of  ground  motion  when  avail- 
able. 

The  variable  frequency  magnitudes  used  for  event 
discrimination  are  based  on  the  envelope  maxima.  As  seen 
in  Figure  3b  (Block  B) , the  traditional  body  and  surface 
wave  magnitude  relationships  are  used  in  the  MARS  program. 

Given  m^  (f)  estimates  at  several  frequencies  within  the 
teleseismic  band  (e.g.,  0.3  to  several  Hertz),  we  can  then 
construct  m^ff^  versus  mb(f2)  plots,  with  f1  <<  f2,  for 
different  sensor-source  region  combinations  and  test  for 
event  discrimination.  Examples  of  this  procedure  as  applied 
to  a large  population  of  Eurasian  events  have  been  described 
in  previous  reports  (Savino,  et  al.,  1975;  Bache,  et  al. , 1976). 
While  m^tf)  versus  Mg(f)  type  discrimination  could  be  easily 
implemented,  as  mentioned  previously  we  are  initially  con- 
centrating on  analysis  of  short-period  P wave  recordings,  or 
m^(f)  type  discrimination. 

2.4  PREVIOUS  DISCRIMINATION  RESULTS 

Most  of  our  previous  experience  with  the  m^ff)  discrim- 
inant has  been  with  LASA  full  array  beam  recordings  of  short- 
period  P waves  from  predominantly  Eurasian  events  (Savino,  et 
al. , 1975;  Bache,  et  al . , 1976).  We  will  briefly  summarize  the 
results  of  these  previous  studies  and  point  out  the  importance 
for  optimum  discrimination  of  determining  the  various  parameters 
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involved  in  the  ( f ) estimates  (e.g.,  the  particular  fre- 

quencies) and  the  application  of  noise  corrections. 

Figure  4 shows  the  discrimination  results  obtained  with 
the  LASA  data.  The  event  data  shown  are  only  two  of  the  fre- 
quency dependent  magnitudes  that  characterize  these  events. 
However,  for  the  range  of  teleseismic  distances  (A  = 60  - 90°) 
and  the  particular  type  of  receiver  (LASA-full  array  beam) , 
these  two  frequency  dependent  magnitudes  at  0.45  and  2.25  Hz 
provide  the  best  separation  of  the  explosion-earthquake  popula- 
tions in  the  m^tf)  parameter  space. 

The  events  in  Figure  4 are  located  in  a variety  of 
geologic  and  tectonic  settings  within  Eurasia:  two  of  the 
presumed  explosions  occurred  in  the  arctic  islands  of  Novaya 
Zemlya;  the  majority  of  the  remaining  explosions  in  the  plat- 
form region  of  eastern  Kazakhstan;  the  earthquakes  originated 
along  the  Alpide  and  western  Pacific  seismic  belts.  Approxi- 
mately 50  percent  of  the  shallow  earthquakes  are  located  along 
the  northern  Japan,  Kuril  and  Kamchatka  arcs.  We  also  point 
out  that  the  Longshot  explosion  which  was  detonated  in  the 
Aleutian  Islands  is  plotted  on  Figure  4 with  coordinates 
(0.45  Hz)  * 5.77  and  (2.25  Hz)  = 4.88.  It  has  been  sug- 
gested in  the  literature  (Von  Seggem  and  Blandford,  1977) 
that  the  locations  of  earthquakes  in  low-Q  (high  attenuation) 
zones,  as  opposed  to  presumed  Eurasian  explosions  which  are 
predominantly  in  high-Q  (shield)  zones,  may  account  for  the 
discrimination  observed  in  Figure  4,  as  well  as  discrimination 
observed  with  other  short  period  discriminants  (Lacoss,  1969). 
While  Q variations  are  important,  and  undoubtedly  contribute 
to  the  scatter  of  the  earthquake  population  in  Figure  4,  the 
results  of  several  different  studies  suggest  that  differences 
in  source  characteristics  are  likely  to  account  for  short- 
period  discrimination.  For  instance,  Barazangi,  et  al.  (1975) 
found  large  lateral  variations  in  P-wave  attenuation  along 
several  of  the  circum-Pacific  arcs.  In  particular,  they  mapped 
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Figure  4.  Variable  frequency  magnitude  estimates  at  fg  = 
0.45  Hz  and  fQ  * 2.25  Hz  based  on  short-period 
P waves  from  Eurasian  events  recorded  at  LASA 
(full  array  beam) . 
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upper  mantle  regions  of  relatively  low  attenuation  (Q  ^ 500) 
beneath  the  shallow  seismicity  along  the  Japan,  Kuril  and 
Kamchatka  arcs.  The  anomalously  low-Q  upper  mantle  zones 
were  located  well  landward  of  the  trench  axes  in  three  regions. 
As  mentioned  previously,  a large  number  of  the  earthquakes  in 
Figure  4 are  located  along  these  arcs  over  the  high  Q zones 
identified  by  Barazangi,  et  al.  (1975).  At  a frequency  of 
2.25  Hz,  the  difference  in  ( 2 . 25 ) for  propagation  through  an 
upper  mantle  with  Q of  500  versus  a Q of  » is  only  0.15  magni- 
tude units,  which  is  less  than  the  separation  of  even  those 
earthquakes  that  plot  closest  to  the  explosion  population.  It 
is  also  hard  to  imagine  that  the  propagation  paths  for  at  least 
a few  of  the  remaining  earthquakes,  that  occurred  in  other 
seismic  zones,  were  not  in  fact  also  along  high-Q  paths. 

An  ideal  experiment  that  would  provide  conclusive 
evidence  for  the  inherent  effectiveness  of  short-period  dis- 
criminants would  be  to  compare  explosions  and  earthquakes  in 
the  same  source  region,  thereby  minimizing  the  contribution 
from  differences  in  propagation  paths.  Two  such  experiments 
were  performed  for  events  at  NTS.  Bakun  and  Johnson  (1970) 
computed  spectral  ratios  based  on  the  short-period  phase 
recorded  at  the  seismograph  station  at  Jamestown,  California 
(JAS) , for  several  NTS  explosions,  afterevents  and  natural 
earthquakes  outside  NTS.  The  distance  range  considered  was 
250  to  500  km.  These  authors  found  that  the  NTS  explosion 
spectra  were  significantly  richer  in  high  frequency  (1.35  to 
2.0  Hz)  than  spectra  for  all  the  natural  earthquakes  (located 
on  both  the  near  and  far  side  of  NTS  with  respect  to  JAS)  and 
possible  afterevents  of  the  NTS  Scotch  explosion.  On  the* 
other  hand,  they  observed  that  the  P-wave  spectra  of  after- 
events of  the  large  NTS  explosions  (Benham,  Boxcar  and  Jorum) 
resemble  the  explosion  spectra. 

At  approximately  -he  same  time,  Basham,  et  al.  (1970) 
reported  the  results  of  spectral  ratio  and  Mg-m^  discrimination 


tests  they  performed  on  Benham  and  several  of  its  afterevents, 
one  of  which  was  included  in  the  experiment  by  Bakun  and 
Johnson  (1970).  Basham,  et  al . (1970)  used  seismograms  written 
at  the  Yellowknife  array  for  their  study.  In  this  case  the 
epicentral  distance  was  approximately  25  degrees.  In  addition 
to  complete  separation  on  an  Mg-mb  plot,  they  noted  that  the 
spectral  ratio  test  produced  excellent  separation  between  the 
larger  Benham  afterevents  and  similar  body-wave  magnitude  NTS 
explosions.  Taken  together,  the  results  of  these  studies  indi- 
cate that  there  are  differences  in  the  short-period  P-wave 
spectra  of  closely  spaced  earthquakes  and  explosions  which  can 
be  exploited  by  various  short-period  discriminants.  The 
apparent  ambiguity  for  the  Benham  afterevents  may  well  be 
related  to  source  mechanism.  While  seen  as  explosion-like  at 
JAS  (A  = 3.5°,  Az  = 285°) , for  a take-off  angle  and  azimuth 
appropriate  for  propagation  to  Yellowknife  (A  = 25°,  Az  = 2°) 
these  events  exhibit  earthquake-like  spectra.  This  potential 
problem  could  be  overcome  by  the  use  of  multi-station  discrimi- 
nation procedures. 

As  the  data  base  for  the  present  discrimination  experi- 
ment increases,  and  when  more  accurate  hypocentral  information 
for  the  events  becomes  available,  we  will  study  our  results  in 
conjunction  with  known  estimated  variations  in  attenuation. 

The  availability  of  data  from  the  numerous  stations  in  this 
experiment  will  hopefully  allow  us  to  differentiate  between 
source  and  path  effects  on  m^ff)  type  discrimination. 

Some  of  the  scatter  in  the  earthquake  and  explosion 
populations  in  Figure  4 is  due  to  background  seismic  noise. 

This,  of  course,  is  relatively  more  important  for  events  of 
small  magnitude  since  the  signal  power  will  be  low  relative 
to  the  noise  power,  and  causes  some  mixing,  or  convergence, 
of  the  two  populations  at  the  very  low  magnitudes.  (The 
smallest  explosions  shown  have  conventional  magnitudes,  m,  , 
of  about  4.) 
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The  LASA  time  series  as  originally  supplied  to  us  did 
not  include  samples  of  the  noise  background  preceding  the  sig- 
nals. In  order  to  study  the  effects  of  applying  noise  correc- 
tions to  the  m^Cf)  data  we  obtained  digital  recordings  of  P 
waves  for  a subset  of  the  events  in  Figure  4 recorded  at  the 
Oyer  subarray  in  Norway.  Each  of  the  time  series  included  at 
least  20  seconds  of  noise  preceding  the  signal.  Figure  5a 
shows  the  pair  of  { f ) estimates  that  provided  the  best  dis- 

crimination for  this  source  region-receiver  combination.  Note 
that  the  magnitude  estimates  in  Figure  5a  are  not  corrected 
for  noise  and,  once  again,  we  observe  convergence  of  the  event 
populations  at  lower  magnitudes.  The  arrows  on  the  smallest 
explosion  data  points  indicate  the  direction  in  which  a noise 
correction  would  move  these  points  in  this  magnitude  plane. 
Figure  5b  shows  the  effect  of  the  noise  correction  applied  to 
all  the  events,  including  earthquakes  as  well  as  explosions. 

(In  this  figure,  only  the  log  of  the  amplitude  is  plotted,  the 
m^  plot  would  be  identical  except  for  a scale  change.)  It  is 
evident  that  the  explosion  population  is  now  well  separated 
from  the  earthquake  population.  The  arrows  on  the  two  small 
explosions,  both  with  conventional  m^  magnitudes  somewhat  less 
than  4.0,  indicate  that  it  is  likely  that  the  unbiased  esti- 
mate of  the  noise  contamination  is  lower  them  that  actually 
present  and  that  the  true  event  points  are  still  lower  in  the 
direction  indicated.  Fortunately,  all  the  data  for  the  pre- 
sent discrimination  experiment  include  adequate  samples  of 
background  noise. 

. The  trends  of  the  earthquake  and  explosion  populations 
in  Figures  4 and  5a  are  approximately  the  same,  even  though 
the  ( f ) planes  are  computed  for  significantly  different  fre- 

quency combinations.  Comparison  of  the  noise  corrected  m^ff) 
estimates  plotted  in  Figure  5b  with  the  m^ff)  planes  in  Figures 
4 and  5a  indicates  that  the  major  effect  of  the  noise  correc- 
tions applied  to  the  Norway  data  has  been  to  rotate  both  the 
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Figure  5a.  Variable  frequency  magnitude  estimates  at  f_ 
0.6  Hz  and  fc  = 5.0  Hz  based  on  short-period 
P waves  from  Eurasian  events  recorded  at  the 
Oyer  subarray  in  Norway.  Arrows  attached  to 
presumed  explosions  indicate  noise  contamina- 
tion of  the  m.  (0.6  Hz)  estimates. 
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earthquake  and  explosion  populations  in  a counterclockwise 
direction.  The  result  is  to  increase  the  slope  of  a line 
separating  the  two  event  populations  in  Figure  5b. 


2.5  PRESENT  DISCRIMINATION  EXPERIMENT 

The  approach  that  we  have  adopted  at  this  early  stage 
in  the  experiment  is  to  build  on  our  previous  discrimination 
results  before  extending  the  analysis  to  include  data  from 
stations  where  we  have  no  experience  with  the  m^ff)  discrimi- 
nant. As  described  in  the  previous  subsection  of  this  report, 
most  of  our  experience  has  been  with  Eurasian  events  recorded 
at  LASA  and  Norway.  In  addition,  a single  presumed  explosion, 
occurring  at  Semipalatinsk  and  recorded  at  the  SRO  station  in 
Albuquerque,  New  Mexico  (ANMO) , has  been  tested  for  discrimina- 
tion with  the  m^(f)  technique  by  Lambert  and  Bache  (1977) . 

A further  criterion  for  the  selection  of  stations  for 
initial  detailed  analysis  is  the  availability  and  quality  of 
event  data,  particularly  for  the  first  10  to  15  events  that 
were  received  earlier  in  the  reporting  period.  Table  2 
summarizes  the  status  of  the  short-period  P wave  data  received 
to  date  on  an  event-station  basis.  We  distinguish  only  three 
possibilities:  + implying  usable  data;  - nonusable  data  and; 

blanks  for  no  data.  Not  included  in  this  table  is  the  fact 
that  many  of  the  event  recordings  at  a particular  station 
have  very  low  signal-to-noise  ratios.  Based  on  Table  2 and 
our  previous  experience,  the  initial  attempts  at  discrimina- 
tion have  concentrated  on  ANMO  and  LAO. 

In  Figure  6 m^f)  estimates  at  2.25  and  0.45  Hz  are 
plotted  for  12  events  recorded  at  ANMO  (closed  circles)  and 
four  events  recorded  at  LAO  (open  circles).  The  m^ff)  esti- 
mates have  been  corrected  for  nominal  instrument  reponses. 

The  dashed  and  solid  curves  represent  the  bounds  of  the 
Eurasian  earthquake  and  explosion  populations,  respectively. 
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Figure  6.  Variable  frequency  magnitude  estimates  at  fc  * 

0.45  Hz  and  2.25  Hz  based  on  short-period  P waves 
recorded  at  ANMO  and  LAO  for  several  of  the 
Eurasian  events  in  the  discrimination  experiment. 
These  estimates  are  compared  to  populations 
defined  in  Figure  4. 
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plotted  in  Figure  4.  Recall  that  the  m^ff)  plotted  in  Figure 
4 were  based  on  the  full  array  beam  recordings  from  LASA, 
whereas  the  open  circles  in  Figure  6 are  magnitude  estimates 
from  the  LAO  subarray  (AO)  recordings.  Finally,  the  solid  tri- 
angle denotes  a presumed  explosion  located  at  Semipalatinsk  and 
recorded  at  ANMO.  This  event  was  the  subject  of  a special 
study  by  Lambert  and  Bache  (1977)  wherein  a multiple  explosion 
scenario,  simulated  with  the  ANMO  time  series,  was  success- 
fully decomposed. 

Considering  the  results  for  ANMO  in  Figure  6,  we  tenta- 
tively identify  events  47  and  49  as  earthquakes  and  events  14, 
19,  20,  21,  22  and  53  as  explosion-like.  Both  the  high  and 
low  frequency  m^  estimates  for  events  16  and  18  are  contami- 
nated by  background  noise  and  will  require  corrections  before 
more  positive  identification  can  be  made.  Note,  however,  that 
at  this  point  these  two  events  are  likely  to  remain  explosion- 
like since  it  was  not  possible  to  make  noise  corrections  for 
the  small  earthquakes  previously  studied  at  LASA  (the  dashed 
population)  with  the  net  effect  that  the  high-frequency  esti- 
mates, in  particular,  were  biased  high.  Had  we  been  able  to 
make  noise  corrections,  similar  to  those  for  Norway  in  Figure 
5b,  we  would  expect  the  low  magnitude  portion  of  the  dashed 
earthquake  population  to  preferentially  move  to  the  left  in 
Figure  6.  Based  on  the  same  argument  we  identify  event  1 as 
explosion  like. 


The  results  for  Station  LAO,  while  more  limited  (see 
Table  2),  agree  with  those  for  ANMO  for  common  events  (19,  47, 
49)  and  indicate  that  event  33  is  explosion-like. 


The  epicentral  distance  range  for  the  events  in  Figure 
6 to  ANMO  and  LAO  is  about  60  to  98  degrees.  Thus,  it  is  of 
interest  to  examine  the  m^ff)  behavior  of  these  events  at 
smaller  epicentral  distances.  Figure  7 shows  magnitude  esti- 
mates for  11  events  recorded  at  the  SRO  station  KAAO  in  Kabul, 
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Figure  7.  Variable  frequency  magnitude  estimates  at  fc  = 
0.6  Hz  and  fc  «*  5.0  Hz  for  a subset  of  the 
Eurasian  events  in  the  discrimination  experiment 
recorded  at  KAAO. 


Afghanistan.  The  distance  range  is  16°  to  65°.  We  selected 
the  nv^(f)  estimates  that  were  least  contaminated  by  noise  and 
span  the  widest  frequency  range  sampled.  Given  these  criteria, 
the  two  narrow-band  filter  center  frequencies  selected  are 
0.6  Hz  and  5 Hz. 

Based  on  the  results  of  ANMO  and  LAO  (Figure  6) , several 
of  the  events  plotted  in  Figure  7 (1,  14,  16,  17,  19,  20,  21, 

22  and  53)  were  tentatively  identified  as  explosion-like;  while 
event  49  was  called  an  earthquake.  Note  that  no  short-period 
P-wave  data  are  available  for  event  47  at  KAAO  (Table  2) . The 
tentative  event  identification  based  on  ANMO  and  LAO  data  is 
supported  by  the  results  from  KAAO  for  all  but  possibly  one 
event,  number  16,  which  is  on  the  border  line  in  Figure  7. 
Obviously,  we  will  have  to  analyze  more  events  at  KAAO  and 
establish  the  boundaries  of  the  earthquake  and  explosion 
populations  before  we  can  make  more  definitive  statements. 

The  next  station  that  we  analyzed  data  from  is  in  Bluff, 
Alaska  (BFAK) . Six  events,  with  epicentral  distances  between 
33°  and  61°  were  available  for  study  and  their  m^ff)  estimates 
at  filter  center  frequencies  of  0.6  Hz  and  4.0  Hz  are  shown  in 
Figure  8.  The  results,  while  restricted  in  number,  are  quite 
striking.  Event  47,  which  was  called  an  earthquake  at  ANMO 
and  LAO,  is  well  separated  from  the  remaining  five  events  which 
appeared  explosion-like  at  ANMO  and  LAO  (Figure  6) . In  addition, 
a noise  correction  for  the  m^  (4.0  Hz)  estimate  for  event  47 
would  further  increase  the  separation.  The  final  point  to  be 
noted  is  the  near  linearity  of  the  explosion-like  events  in 
Figure  8.  We  will  be  particularly  interested  in  seeing  how 
the  populations  fill  in  at  this  station  as  more  events  are 
analyzed. 


37 


m^  [4 . 0 Hz] 


Figure  8 


Variable  frequency  magnitude  estimates  at  f 
0.6  Hz  and  fc  = 4.0  Hz  for  a subset  of  the 
Eurasian  events  recorded  at  BFAK. 


2.6 


SUMMARY 


During  this  reporting  period,  we  received  seismic  wave- 
forms for  nineteen  Eurasian  events  and  initiated  an  event 
identification  study  based  on  variable  frequency  magnitude 
estimates.  Table  3 summarizes  the  preliminary  identifications 
made  for  thirteen  of  the  larger  events  recorded  at  four  dif- 
ferent stations.  At  this  point  in  the  experiment  we  are  unable 
to  assign  probabilities  to  the  event  identifications  in  Table 
3.  To  do  so  will  require  the  addition  of  more  events  so  that 
the  extent  of  the  earthquake  and  explosion  populations  can  be 
better  defined  at  each  of  the  stations. 

Our  plans  for  the  future  are  the  following: 

1.  Define  the  event  populations  at  all  the  stations 
contributing  short-period  P-wave  seismograms. 

2.  Apply  noise  corrections  to  the  m^(f)  estimates. 

3.  Test  automated  single-  and  multi-station 
discrimination  techniques  once  the  data  sets 
are  adequate. 
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TABLE  3 

PRELIMINARY  EVENT  IDENTIFICATION 
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Earthquakes  Explosion-Like  Events 


(Event 

Number) 

(Event 

Number) 

ANMO 

LAO 

KAAO 

EURO 

ANMO 

LAO 

KAAO 

EURO 

47 

47 

47 

1 

1 

49 

49 

49 

14 

14 

16 

16 

17  17 


18 

19  19  19 

20  20 

21  21  21 

22  22 

33  33 

53  53  53 
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III.  ANALYTIC  CONTINUATION  OF  THE  ELASTIC  FIELD  FROM 
A COMPLEX  SOURCE  IN  A HALFSPACE 

An  important  problem  in  nuclear  explosion  and  earthquake 
seismology  is  the  development  of  a capability  for  linking 
detailed  numerical  source  calculations  with  efficient  analytical 
techniques  for  propagating  elastic  waves  in  layered  media. 

If  the  source  is  in  a homogeneous  wholespace,  techniques 
employing  an  expansion  of  the  outgoing  waves  in  spherical 
harmonics  work  very  well  and  have  been  employed  to  study 
earthquakes  (Bache,  et  al. , 1976)  and  complex  nuclear  explo- 
sions (Cherry,  et  al. , 1975,  1976;  Bache  and  Masso,  1978). 
However,  these  techniques  are  not  rigorously  valid  when  a 
material  boundary  is  present  in  the  source  region.  Bache, 
et  al.  (1977)  introduced  some  ad  hoc  approximations  and 
applied  the  whole  space  technique  to  study  near-surface 
explosions  in  a halfspace.  Even  though  this  worked  fairly 
well  for  the  case  studied,  it  has  some  serious  drawbacks 
and  better  methods  are  needed. 

In  this  section,  we  present  some  mathematical  results 
outlining  a method  for  analytically  continuing  the  output  of 
a finite  difference  simulation  of  a source  in  a layered  half- 
space. The  source  calculations  may  include  arbitrary  non- 
linearity in  a portion  of  the  grid.  It  is  also  necessary 
that  the  waves  be  computed  out  into  the  region  where  the  mate- 
rial response  is  linearly  elastic. 

The  particular  case  discussed  here  is  the  computation 
of  the  body  waves  at  an  arbitrary  range  in  a layered  earth 
model  from  an  axisymmetric  source  in  a half space.  Following 
a similar  procedure,  the  surface  waves  can  also  be  computed, 
though  this  is  not  worked  out  here.  The  method  used  may  be 
generalized  to  other  geometries. 
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The  mathematical  basis  of  our  approach  is  the  elasto- 
dynamic  integral  representation  given  by  Burridge  and  Knopoff 
(1964).  This  is  a rigorous  expression  for  the  displacement  in 
an  elastic  medium  in  terms  of  a surface  integral  over  a (non- 
physical) boundary  E which  entirely  encloses  the  source  volume 
(including  all  inelastic  regions) . Two  elements  are  required 
to  form  this  integral  solution.  First,  we  must  have  the  time 
histories  of  the  displacement  vector  U (xQ,t)  and  the  stress 

vector  r (x  . t)  on  E.  These  are  obtained  from  the  finite  dif- 
— o 

ference  source  calculation.  Second,  the  Green's  tensor 

g (x,  x ; t-t  ) and  its  spatial  derivatives  (for  the  appropriate 
® ® 

elastic  medium)  must  be  evaluated  for  the  prospective  receiver 
point  x,  for  all  ^ on  I.  The  elastic  medium  used  to  compute 
U (xQ,t)  and  x (xQ,t)  must  be  the  same  as  that  used  to  compute 
g (x,  3^;  t-tQ). 

The  general  form  of  the  integral  representation  (assuming 
isotropy)  is 


U (x,t)  = J dt  Jds  [g (x,  Xq;  t- 


V • i <*o'  V 


- g Z t-tQ)  : M (Xo»t0)] 


(3.1) 


where 


M(x  ,t  ) = XI  n • U(x  ,t  ) + u [n  U(x  ,t  ) + U(jc  ,t  )n]. 


x and  t are  receiver  position  and  time, 

X and  u are  elastic  constants, 

n is  the  normal  to  E,  directed  into  the  source  volume, 
U is  the  displacement  vector, 
jr  is  the  stress  vector  on  E, 
g is  the  Green's  tensor  solution, 

I is  the  identity  tensor. 
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This  formulation  is  entirely  general.  However,  we  will 
now  specialize  to  the  case  of  an  axisymmetric  source  and  assume 
that  the  geologic  structure  is  represented  by  some  plane- 
layered earth  model. 

We  assume  that  the  elastic  medium  occupies  the  halfspace 
z > o.  We  introduce  the  cylindrical  coordinates  r,  9,  z,  and 
assume  that  the  elastic  properties  of  the  medium  and  the  r and 
z components  of  displacement  are  independent  of  9 and  that  the 
0 component  of  displacement  is  zero.  Finally,  we  specify  the 
surface  E to  be  a cylinder  of  radius  a extending  from  z * 0 
to  z = b.  The  geometry  is  shown  in  Figure  9.  In  this  special 
case.  Equation  (3.1)  can  be  written  in  terms  of  cylindrical 
components  as 


(T  (r , z , t)  = 


a CGr  t,z  *rr 


G3*  „ a ] 

z t,z  rz 


- tGr*t,r  arz  * Gz*t,r  °zz) 


where  G?  is  the  azimuthally  averaged  Green's  tensor  component 
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2 TT 

Gi*r0'r'z0'z't“t0^  = j gi^r0,r,z0,z  <^“^0 ' t— fc0  ^ d<i>  (3*3^ 

0 

and  the  operators  * , and  * , are  defined  by 

tl  Z u IT 


G3  * f 
i t,z 


/ dt°  1 


d2n  (a,r,z_ ,z,t-t.)  f (a,zrt,t-) 


•o  ri 


O'  O' 


® a 

Gi*t,rf  = A dt0  /"  <^0  jrQG?  (rQ  ,r,b,z,t-tQ)  f (r0,b,tQ)j 

^ M J r\ 


(3.4) 


Equation  (3.2),  while  lengthy,  is  a sura  of  similar  terms, 
Each  term  has  the  form 


Gi(r'ro,zo'z't“to)  *t'Z  vi(ro'zo,to) 


where  v^  is  a displacement  or  traction  component  monitored  on 
Z.  Thus,  the  procedure  is  to  solve  for  the  Green's  function 
for  a point  on  E,  carry  out  the  integrals  (Equation  3.4)  with 
the  stress  and  displacement  time  histories  at  that  point  and 
then  sum  the  integrals  according  to  Equation  (3.2). 

The  form  of  the  Green ' s functions  depends  on  the  range 
and  type  of  seismic  wave  to  be  propagated.  In  the  following, 
we  present  the  Green's  functions  that  are  compatible  with  far- 
field  surface  wave  methods  and  with  near-field  and  far-field 
body  wave  methods. 

We  begin  with  the  integral  solutions  for  the  displace- 
ment due  to  a point  force  embedded  in  a homogeneous , linearly 
elastic  space: 

GO 

g3(r,z,y)  - ■ -1  2 f dk  jk  [A* . (k,u)  +A*  (fc,u)]| 

47rpui  J.  3 ^ 
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(3.5) 


t 


where 


Aa  . = — 
ij  v 


-V  z 

a 


r,.2 


KJ0(kr) 


-kv  J. (kr) 
a 1 


-kvoJl(kr) 

-vaJ0  <kr>  . 


A3. 

i: 


-v 


r-VgJQ(kr)  kv3J1C<r)' 


Lkv3JL(kr)  k JQ (kr) 


(3.6) 


and 


Y 


k2  - - 

/ Y 


2\  1/2 


The  azimuthally  averaged  Green's  tensor,  Equation  (3.3),  is 
required.  As  is  shown  in  Figure  10,  we  introduce  the  dimen- 
sions a,  the  cylinder  radius  and  rT,  the  horizontal  distance 
from  the  center  of  the  cylinder  to  the  observer,  while  r is 
the  horizontal  distance  from  a point  on  the  cylinder  to  the 
observer.  The  only  functions  in  Equation  (3.5)  which  depend 
on  r are  the  Bessel  functions.  Now,  from  Watson  (1944), 


Jv (kr) e 


ViO 


£ 


m=- 


< . % _ m % « ojno 

Ja  jCa  e 


u 


c 


and,  taking  the  azimuthal  average,  we  have 


2t 


Jv (kr) d$ 


2 W**’  Vka)  / 

m»- « 0 


J2v(krT)  Jv(ka) . 
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Figure  10.  The  monitoring  cylinder  I is  shown  from 
directly  above  the  free  surface.  The 
geometry  for  computing  the  azimuthal 
average  of  the  Green's  function  is  also 
shown. 
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The  azimuthally  integrated  Green's  tensor  is  then 


G?  (w) 


where 


CO 

if 


dk  k 


3°\  + 3?  . ] j ' 

ij  iJJ  )- 


(3.7) 


r,  2 


-v  z 
Ba  = 6 -.2- 

Bij  6 


k‘JQ(krT)  JQ(ka)  -kvaJ2 (krT) J± (ka) 
[-kvaJ2  (krT)  Jx  C«)  -v|J()  (krT)  J0  (ka) 


ana 


r ..2 


3ij  * - 


e 


“V 


3 


-v3JQ (krT) JQ (ka)  kvgJ2 (krT)  (ka) 


2. 


[_kvg  J2  (krT)  JL  (ka)  k JQ  (krT)  JQ  (ka)  J 


Equation  (3.7)  may  be  combined  with  layer  matrices  to  develop 
the  theory  for  surface  waves. 

To  compute  body  waves,  we  have  chosen  to  use  generalized 
ray  methods  and  the  Cagniard-de  Hoop  inversion  technique.  To 
cast  the  solutions  in  the  form  compatible  with  these  methods, 
we  return  to  Equation  (3.5).  We  make  the  substitutions  oj  = -is 
and  k = -isp.  This  separates  the  s and  p dependence  and  takes 
us  from  Fourier  to  Laplace  transform  space.  The  solutions  for 
the  point  forces  are  now 

i« 


g? (r,z,s) 


2pt 


la 


/ 


do  c 


-ST]  z 


CTj (p,s) 


-ST^z  3 

+ « " C?  j (p , s) 
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(3.8) 


» 


where 


C“  (p,s)  .i- 

CL 


rp2KQ (spr) 


~£pnaKl ^sPr) 


_-£napK1(spr)  naXQ(spr)  J 


C^fp.s,  - J- 


HgKg  (spr)  eprigK1  (spr) 
L £PHgKi (spr)  p2XQ(spr)  J 


ana 


C 


S = 


■fr- 


1/2 


+1  for  upgoing  wave. 


-1  for  downgoing  wave. 


It  can  be  shown  (e.g.,  Cistemas,  Betancourt  and 
Leiva,  1973)  that  the  response  of  a layered  medium  with  L 
layers  with  thickness  h^  and  elastic  properties  a^,  3^/ 
is  given  by 


gi(r,z,s) 


s 

2 PTT* 


ao  j^ao 

E -/ 


dp  p 


n=l 


+C?j(p,s)  T® (p) exp (-st®) 


Cij(p,s)  Tn(p)  exp(-st~) 


(3.9) 


Each  element  of  the  sum  in  Equation  (3.9)  represents  a "gener- 
alized ray"  prescribed  by  the  layers  a ray  traverses  and  its 
mode  (p  or  s)  in  each  layer.  A ray  has  a vertical  phase  shift 
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*n  * (nTh)s  *iLlhlnyi 
SL- 1 

(the  subscript  s denotes  properties  in  the  source  layer) 
and  amplitude  factor  TY,  the  product  of  reflection  and  trans- 
mission coefficients  due  to  the  impedance  contrasts  between 
the  layers.  Although  the  sum  in  Equation  (3.9)  is  infinite 
for  problems  related  to  the  work  proposed  here,  the  solution 
is  dominated  by  a few  rays. 

Before  we  can  use  Equations  (3.8)  and  (3.2),  we  must 
perform  an  inverse  Laplace  transform  and  compute  the  azimuthal 
average.  First  transforming,  we  obtain  the  time-domain  Green's 
tensor  for  the  n^*1  generalized  ray.  This  is 


g^  (r , z , t)= 
n 


^dTpC“j(P,t) 
' fcon 


T^(P)  [ ( t— t ) (t— x+2pr) ] 1/2 

n q c 


+ Im  f dtp  C®.(p,t)  T*(p)  [(t— T) (t— x+2pr) ] 1/2 

in 


(3.10) 


where  the  path  of  integration  is  defined  by 


Im  t =*  Im 


N 

?r  + S V*  + hsns 
£=1 


0, 


(3.11) 


and  t*  is  the  ray  arrival  time.  The  matrices  C in  the  time 
on 

domain  are  given  by 


I 
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m fa- 


I 


5 (P/t) 


t-T+or 

pr 


In  the  near-field  (where  the  periods  of  interest  and 
the  travel  time  are  of  the  same  order) , the  integrals  in 
Equation  (3.10)  must  be  done  numerically . The  results  for 
different  azimuths  (<j>)  must  then  be  summed  to  compute  the 
azimuthal  averages. 

In  the  far-field,  the  computation  of  (s)  can  be 
simplified  by  invoking  the  approximation  t-T  <<  pr.  In 
this  case  £(p,t)  a 1,  and  Equation  (3.10)  can  be  written  as 


gj  (t) 


+ la  | /p  C 

where  the  functions  of  p are  evaluated  along  the  Carniard 
contour  (Equation  3.11).  The  convolution  may  now  be  done 
after  all  the  rays  are  summed,  which  greatly  reduces  the 


««*•*>  HI  **= 


(3.12) 


computing  time. 

The  far-field  regime  allows  another  time  saving  approx- 
imation in  addition  to  the  one  just  described.  At  distances 
large  compared  to  the  cylinder  dimensions,  the  amplitude  and 


frequency  content  of  Green's  functions  due  to  rays  coming  from 
» various  parts  of  Z are  nearly  equal.  The  greatest  variation 

is  in  their  relative  arrival  times.  Thus,  it  is  necessary  to 
compute  only  one  or  perhaps  a few  Green's  functions  rather 
than  one  for  each  point  on  I.  This  feature  makes  the  azimuthal 
• average  quite  easy  to  compute. 

The  theory  presented  above  is  now  being  programmed  and 
tested.  With  the  program  we  will  be  able  to  continue  the  out- 
put of  axisymmetric  finite  difference  calculations  to  any  range 
of  interest  in  plane-layered  earth  models.  The  analytic  con- 
tinuation is  done  entirely  in  the  time  domain,  obviating  the 
need  to  compute  Fourier  transforms  which  can  sometimes  be 
rather  difficult.  A similar  theory  can  be  developed  for  com- 
puting the  surface  waves  excited  by  the  finite  difference  out- 
put. Having  the  results  presented  here,  the  development  of 
this  theory  is  rather  straightforward  and  will  be  done  during 
this  quarter. 
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IV.  THEORETICAL  COMPUTATION  OF  Lg 


In  this  section  we  summarize  the  results  of  a brief 
study  of  the  excitation  of  the  regional  phase  Lg.  The  study 
was  conducted  at  VSC  request  and  the  results  were  presented 
orally  at  an  ARPA  Program  review  covering  seismic  detection 
and  discrimination  at  regional  distances  and  held  in  Dallas 
on  11  and  12  January  1978. 

Knopoff,  Schwab,  Panza  and  several  others  at  UCLA  and 
elsewhere  have  demonstrated  that  waves  having  the  observed 
characteristics  of  Lg  are  present  in  plane-layered  models 
(e.g.,  Panza,  Schwab  and  Knopoff,  1972;  Knopoff,  Schwab  and 
Kausel,  1973;  Panza  and  Calcagnile,  1975) . Higher  modes  for 
both  Rayleigh  and  Love  waves  in  continental  structures  have 
stationary  phases  at  about  3.5  km/sec. 

To  illustrate  the  theoretical  behavior  of  Lg,  synthetic 
seismograms  were  computed  for  a double-couple  source  in  an 
earth  model  meant  to  represent  an  average  continental  structure. 
The  earth  model  was  provided  by  VSC  (Michael  J.  Shore,  per- 
sonal communication) . The  shear  wave  velocity-depth  profile  for 
this  model  is  plotted  in  Figure  11.  Also  plotted  in  the  figure 
is  the  low  velocity  zone  (LVZ)  model  for  which  the  Lg  phase 
was  studied  extensively  by  Knopoff,  Schwab  and  Kausel  (1973) , 
Knopoff,  et  al.,  (1974)  and  Panza  and  Calcagnile  (1975).  The 
most  important  difference  between  the  two  is  that  the  VSC  model 
has  no  sedimentary  layers  at  the  top.  Instead,  the  top  12 
kilometers  are  uniform  with  6 * 3.58  km/sec.  Otherwise,  the 
two  models  are  sufficiently  similar  for  the  trends  reported  by 
Knopoff  and  coauthors  to  remain  valid.  Also  provided  by  VSC 
was  a Q model  which  takes  a simple  form.  The  Q is  1500  to  a 
depth  of  300  km  and  increases  at  greater  depth.  For  our  pur- 
pose this  is  essentially  a constant  Q model. 


The  Love  wave  phase  and  group  velocities  and  source-depth 
independent  amplification  A^  (Harkrider,  1964)  for  the  VSC  struc- 
ture of  Figure  11  are  plotted  in  Figure  12  for  periods  down  to 
one  second.  There  is  also  a depth  and  frequency  dependent  term 
which  is  not  plotted  here.  In  Figure  12  the  group  velocity 
stationary  phase  at  about  3.5  km/sec  is  clearly  apparent. 
Associated  with  the  sudden  drop  in  the  group  velocity  from  4.5 
to  3.5  km/sec  is  a sudden  increase  in  the  A^  Thus,  the  higher 
group  velocity  waves  are  much  less  strongly  excited. 

Similar  plots  are  shown  for  Rayleigh  waves  in  Figure  13. 
These  plots  point  out  a difficulty  with  the  structure  used  for 
the  calculations.  At  short  periods  the  fundamental  mode  has  a 
group  velocity  that  is  not  much  smaller  than  3.5  km/sec.  Adding 
slower  sedimentary  layers  near  the  surface  would  cause  the  funda- 
mental mode  to  drop  to  substantially  lower  velocities  at  short 
periods. 

Some  typical  theoretical  seismogram  calculations  for 
transverse  Lg  are  shown  in  Figures  14  and  15.  The  modes  com- 
puted are  those  for  which  the  phase  and  group  velocities  and 
amplification  are  shown  in  Figure  12.  The  source  is  a strike- 
slip  double-couple  and  the  azimuth  is  22.5°  from  the  strike. 

The  ground  motion  is  filtered  by  a WWSSN  long  period  seismom- 
eter which  peaks  at  15  seconds.  On  each  seismogram  is  given 
the  peak-to-peak  amplitude  in  microns  at  15  seconds. 

In  Figure  14  we  compare  the  Lg  at  two  distances  and  two 
source  depths.  We  notice  that  at  these  shallow  depths  the 
seismogram  is  dominated  by  the  fundamental  mode.  The  ampli- 
tudes of  the  higher  modes  decrease  rapidly  with  increasing  mode 
number.  The  rate  of  attenuation  of  the  peak  amplitude  for  each 
mode  between  these  two  distances  is  shown  on  the  figure.  The 
distance  attenuation  is  generally  more  rapid  for  the  higher 
modes,  though  this  is  not  always  the  case  for  these  examples. 
Unfortunately,  the  final  summed  seismograms  do  not  look  very 
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Figure  12.  The  Love  phase  and  group  velocities  and  depth-independent  amplification  are 
plotted  for  the  VSC  provided  structure  of  Figure  13.  The  fundamental  (0) 
and  first  three  higher  modes  are  plotted. 
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Figure  15. 
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Transverse  Lg  component  seismograms  are  compared 
for  three  source  depths  at  a range  of  1000  km. 
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much  like  observed  data.  Much  of  this  can  be  attributed  to 
the  overly  simple  crustal  model  used. 

In  Figure  15  we  compare  the  transverse  Lg  seismograms 
at  1000  km  for  three  depths.  At  33  km  depth  the  fundamental 
mode  no  longer  is  the  largest  and  the  four  modes  computed  are 
roughly  the  same  size.  The  summed  seismogram  is  a bit  more 
complex,  but  still  bears  little  resemblance  to  observations 
of  the  surface  waves  from  earthquakes. 

The  actual  source  used  for  the  seismograms  of  Figures 
14  and  15  was  the  double-couple  contribution  to  a composite 
source  made  up  of  an  explosion  plus  a tectonic  release  com- 
ponent. The  ratio  between  the  two  components  is  that  derived 
by  Toksoz  and  Kehrer  (1972)  for  the  PILEDRIVER  event;  that  is, 
F = 3,  where  the  F factor  represents  the  ratio  between  the 
explosion  and  double-couple  portions  of  the  field.  With  F 
this  large,  the  double-couple  radiation  dominates  at  most 
azimuths. 

In  Figure  16  we  show  the  Rayleigh  waves  from  the  com- 
posite source.  Again,  the  azimuth  is  22.5°  from  the  strike  of 
the  double-couple.  They  are  compared  to  the  Rayleigh  waves 
from  the  explosion  (spherically  symmetric)  portion  of  the 
field  alone.  Addition  of  the  double-couple  makes  some  dif- 
ference in  the  waveforms.  However,  in  either  case  the  higher 
modes  make  almost  no  contribution. 

The  examples  in  Figures  14  through  16  are  not  particu- 
larly interesting  because  of  the  overly  simple  crustal  model 
used  in  the  calculations.  Therefore,  the  scaling  properties 
of  the  Lg  phase  on  these  theoretical  records  are  not  very  use- 
ful. However,  they  do  demonstrate  the  capability  to  compute 
synthetic  seismograms  at  this  range.  This  capability  could 
be  exercised  with  more  realistic  earth  models  to  determine 
the  theoretical  dependence  of  Lg  on  the  source  and  path 
parameters . 
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V.  ANALYSIS  OF  SURFACE  WAVES  FROM  THE  FRENCH 
TEST  SITE  IN  THE  SAHARA 


5.1  INTRODUCTION 

Surface  wave  recordings  from  stations  in  and  near 
Africa  were  used  to  study  the  French  Sahara  underground  event 
SAPHIRE  (27  February  1965) . The  objective  of  this  study  is  to 
better  understand  the  factors  that  control  the  amplitude  and 
waveform  of  explosion-generated  surface  waves  at  regional  and 
teleseismic  distances.  A thorough  understanding  of  these  fac- 
tors is  required  if  surface  wave  signatures  are  to  be  a 
reliable  source  of  information  about  seismic  source  parameters. 

The  approach  taken  in  this  study  is  the  same  as  in  an 
earlier  study  under  this  contract  (Bache,  Rodi  and  Harkrider, 
1978)  in  which  we  modeled  Rayleigh  waves  from  NTS  explosions. 

To  model  the  effects  of  earth  structure,  we  invert  the  ob- 
served dispersion  for  each  path  to  obtain  a model  of  the 
average  crust  and  upper  mantle  structure  along  the  path. 

These  path  models  and  a model  of  the  structure  near  the 
source  should  largely  account  for  the  waveforms  of  the  ob- 
served surface  waves  if  the  theory  we  use  is  appropriate.  To 
test  this,  we  compare  the  observed  waveforms  to  synthetic 
seismograms  generated  from  our  earth  models  and  a trial  model 
of  the  source.  If  the  waveforms  are  matched  successfully, 
we  can  then  interpret  the  surface  wave  amplitudes  at  the 
stations  in  terms  of  the  source  function. 

This  approach  and  its  theoretical  basis  are  described 
in  more  detail  in  the  next  subsection.  Subsequent  subsections 
present  the  surface  wave  data,  inversion  results  and  seismo- 
gram synthesis  -results  for  the  SAPHIRE  event. 

5 . 2 MODELING  PROCEDURE 

The  time-domain  amplitude  and  waveform  of  surface 
waves  depend  on  the  characteristics  of  both  the  source  and 
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propagation  path.  To  model  this  dependence  we  use  a modifica- 
tion of  Harkrider's  (1964,  1970)  theory  of  surface  wave  genera- 
tion by  point  sources  in  plane-layered  media.  The  modification 
accounts,  in  an  approximate  way,  for  the  transmission  across 
the  boundary  between  the  source  region  and  the  rest  of  the  path 
to  the  receiver.  Thus,  the  earth  may  be  approximated  by  two 
one-dimensional  models,  the  first  representing  the  structure 
near  the  source  and  the  second  representing  the  average  struc- 
ture along  most  of  the  path. 

According  to  the  theory,  the  major  contributions  to 
surface  wave  spectra  can  be  broken  down  as  follows  (see 
Harkrider,  1964,  1970  and  Bache,  Rodi  and  Harkrider,  1978  for 
notation  and  further  details) : 


Amplitude  Spectrum 


Al.  Source-independent  amplitude  response  of  path 

A^  for  Love 

waves)  . 

-Y2r 

A 2.  Dissipative  attenuation  along  path  (e  ). 


structure  (A^  f°r  Rayleigh  waves. 


A3.  Source-independent  amplitude  response  of  source- 
region  structure  (AR1  and  A^)  . 

A4.  Source-depth-dependent  modal  excitation  (K-.  for 
explosive  source) . 


A5.  Radiation  pattern  of  source  (depends  on  source 
size,  time  history  and  geometry) . 


Phase  Spectrum 


PI.  Dispersion  along  path  (ur/c2) . 

P2.  Initial  phase  of  source  (depends  on  same  factors 
as  A5) . 


These  phenomena  are  most  easily  understood  in  the  frequency 
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domain  where  they  are  nearly  separable.  To  understand  them 
in  the  time  domain,  we  must  rely  on  Fourier  transform  theory. 

The  biggest  task  in  our  modeling  procedure  is  the 
estimation  of  the  average  crust  and  upper  mantle  structure  for 
each  source-receiver  path.  This  is  accomplished  by  inverting 
the  observed  dispersion  for  each  path  using  other  data  such  as 
refraction  results  as  constraints.  With  certain  assumptions 
about  the  initial  phase  of  the  source,  path  dispersion  can  be 
determined  nearly  independently  of  other  contributions  to  the 
surface  wave  spectra.  Specifically,  we  assume  small  phase 
and  group  delay  at  the  source.  Rayleigh  and  Love  wave  group 
velocities  are  then  obtained  by  narrow-band  filtering  and 
Hilbert  transform  techniques.  To  determine  Rayleigh  wave 
phase  velocities,  we  assume  the  spectra  are  dominated  by  the 
explosive  component  of  the  source  which  has  a very  small 
initial  phase.  However,  the  phase  spectrum  is  indeterminate 
by  a multiple  of  2ir  and  there  are  a multiplicity  of  phase 
velocity  curves  for  each  path.  To  determine  Love  wave  phase 
velocities  requires  a more  complete  knowledge  of  the  source 
radiation  pattern. 

We  invert  the  phase  and  group  velocity  data  for  each 
path  by  formal  linear  inversion  techniques  to  determine  models 
of  seismic  velocity  and  density.  A first  test  of  the  appli- 
cability of  our  theory  is  whether  reasonable  plane-layered 
models  can  be  found  that  fit  the  observed  dispersion  and 
whether  the  models  for  different  paths  describe  a consistent 
picture  of  the  entire  region  of  study. 

A source-region  model  is  derived  from  one  or  more  of 
the  path  models,  modified  as  necessary  to  agree  with  any 
available  information  about  the  near-surface  structure  of 
the  source  region.  As  a result,  the  models  we  invert  from 
the  observed  dispersion  (PI)  are  expected  to  account  for  spec- 
tral contributions  Al,  A3  and  A4  as  well.  (The  depth  of  the 


source  is  assumed  known.)  A Q or  y model  is  determined  from 
global  or  regional  (if  available)  attenuation  data  and  is  used 
to  account  for  A2.  These  contributions  (all  but  A5  and  P2) 
determine  the  waveform  of  the  surface  waves.  Using  our  path 
and  source-region  models  plus  a trial  model  for  the  source, 
we  synthesize  seismograms  for  comparison  with  the  observa- 
tions. The  agreement  in  waveform  is  a test  of  our  modeling 
procedures  and  the  theory  they  are  based  on.  Furthermore,  if 
our  modeling  is  adequate,  the  observed  and  synthetic  amplitude 
spectra  at  each  station  should  differ  at  long  periods  only  by 
a constant  factor.  The  values  of  this  factor  at  the  various 
stations  can  be  used  to  infer  the  difference  between  the  true 
source  parameters  and  the  parameters  of  the  trial  source 
model . 


5.3  SURFACE  WAVE  OBSERVATIONS 

Figure  17  shows  the  location  of  the  Hoggar  (HOG)  test 
site  in  Algeria,  together  with  the  paths  to  eleven  WWSSN 
stations  in  and  near  Africa.  Reported  here  are  the  results 
from  modeling  the  SAPHIRE  surface  waves  recorded  at  four 
of  these  stations  (HLW,  SHI,  AAE  and  NAI)  plus  another  not 
plotted  (X) . Table  4 gives  a description  of  the  five  paths 
studied.  In  Figure  18  are  plotted  the  seven  seismograms 
analyzed  for  this  report.  These  include  the  vertical  component 
from  each  station  and  the  transverse  component  from  HLW  and 
AAE.  All  of  these  data  were  hand-digitized  from  film-chip 
reproductions.  The  transverse  component  data  were  obtained  by 
rotating  the  horizontal  recordings  with  respect  to  the  great 
circle  path  from  source  to  receiver. 

The  dispersed  surface  wave  is  easily  identifiable  on 
most  of  the  seismograms  in  Figure  18,  but  some  are  clearly 
contaminated  by  noise  or  other  arrivals.  In  particular  the 
AAE  transverse  component  is  very  complex.  Also,  the  Rayleigh 
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Figure  17.  The  location  of  the  Hoggar  test  site  and  the 
paths  to  eleven  WWSSN  stations  are  shown  in 
a polar  projection  centered  on  the  test  site. 
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Figure  18.  The  SAPHIRE  seismograms  analyzed  in  this  study 
are  shown,  with  tQ  denoting  the  start  time  (in 
seconds)  of  each  record  relative  to  the  origin 
time. 
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wave  has  obviously  leaked  onto  the  HLW  transverse  component, 
suggesting  that  the  direction  of  propagation  for  the  HOG-HLW 
path  is  not  great-circle.  Another  example  of  contamination 
is  the  very  long-period  trail  behind  the  peak  arrival  on  the 
X vertical  component. 

Group  velocities  were  determined  from  all  seven  seismo- 
grams with  narrow-band  Gaussian  filtering.  For  the  Rayleigh 
waves,  clear  fundamental-mode  group  arrivals  could  be  obtained 
at  periods  as  low  as  5 and  as  high  as  35  s,  but  mainly  in  the 
7 to  25  s range.  Love  wave  group  velocities  for  HLW  and  AAE 
were  obtained  from  12  to  25  s.  Figure  19  compares  the  group 
velocities  determined  for  the  five  paths  and  used  in  the  inver- 
sion for  the  path  structures. 

The  determination  of  phase  velocities  presented  much 
more  difficulty,  but  was  pursued  because  they  add  substantially 
to  the  resolving  power  for  earth  structure.  Phase  velocity  (c) 
is  derived  directly  from  the  surface  wave  phase  spectrum  (<M 
through  the  relation 

= *s(w)  - 4>  (o>)  + n2ir,  (5.1) 

where  r is  the  source-receiver  distance  and  <j>  is  the  initial 

s 

phase  of  the  source,  including  radiation  pattern  effects.  For 
Rayleigh  waves  we  assumed  <p  = -3u/4,  a close  approximation  to 
the  initial  phase  of  the  vertical  Rayleigh  wave  displacement 
due  to  an  explosive  source.  In  the  case  of  Love  waves,  we  did 
not  feel  confident  in  assuming  <J>  , so  no  Love  wave  phase  velo- 
cities  were  used  for  the  inversions  reported  in  this  report. 

A difficulty  with  deriving  phase  velocity  from  Equation 
(5.1)  is  the  possible  contamination  of  the  phase  spectrum  by 
noise,  body  waves  or  higher  modes.  This  problem  was  handled 
by  time-domain  windowing  before  obtaining  the  Fourier  spectrum 
in  order  to  insure  that  the  spectrum  at  the  periods  of  interest 
(5  to  30  s)  was  dominated  by  fundamental-mode  Rayleigh  wave. 
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Figure  19.  The  Rayleigh-wave  group  velocities  (Ur) 
determined  for  five  paths  from  Hoggar 
test  site  are  plotted  as  open  symbols. 
Love-wave  group  velocities  (Ul)  are 
plotted  as  solid  symbols. 
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We  checked  this  by  comparing  the  group  velocities  implied  by 
our  phase  velocity  estimates  with  the  narrow-band  group 
velocities.  (This  can  be  done  independently  of  the  value  of 
n in  Equation  (5.1)).  Phase  velocities  at  periods  for  which 
the  agreement  was  poor  or  for  which  the  spectral  amplitude  was 
small  were  rejected. 

A second  problem  with  determining  phase  velocities  is 
the  fact  that  the  integer  n in  Equation  (5..1)  is  indeterminate 
from  the  seismograms  alone.  Other  information  is  needed  to 
ascertain  which  of  the  possible  phase  velocity  curves  implied 
by  Equation  (5.1)  is  most  nearly  correct.  For  the  long  paths 
and  short  periods  available  for  this  study,  the  possible  curves 
do  not  differ  a great  deal.  However,  we  were  able  to  choose  n 
with  some  confidence  for  the  HOG-X  path.  Implausible  mantle 
shear  velocities  (>  5 km/s  at  70  km  depth)  resulted  from 
inverting  the  HOG-X  phase-velocity  curve  adjacent  to  the  one 
chosen  (n  differing  by  1).  For  the  other  paths,  we  chose  n 
such  that  the  phase  velocity  curves  for  all  the  paths  were 
most  nearly  aligned  at  the  longest  available  periods. 

The  final  phase  velocities  determined  for  the  five 
paths  are  plotted  in  Figure  20.  Table  5 gives  the  approximate 
shifts  to  these  data  that  result  from  altering  n by  + 1 from 
the  values  chosen.  Our  final  phase  velocity  data  are  con- 
sistent with  long  period  surface  wave  observations  reported 
from  previous  studies  of  east  African  paths.  Knopoff  and 
Schlue  (1972)  found  c(20s)  =3.6  km/s  for  a path  from  AAE  to 
NAI.  Gumper  and  Pomeroy  (1970)  found  c(30s)  = 3.92  for  a 
north-south  path  from  HLW  to  Lwiro,  Republic  of  the  Congo 
(west  of  NAI) . 
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VELOCITY  (km/s) 


Figure  20.  The  observed  Rayleigh-wave  phase  velocities 
(CR)  are  shown  for  five  paths  from  Hoggar 
test  site. 


TABLE  5 


APPROXIMATE  PHASE  VELOCITY  UNCERTAINTIES 
DUE  TO  2tt  AMBIGUITY  IN  PHASE 


Phase  Velocity  Uncertainty  (km/s) 


Station  Period  = 
HLW 
SHI 
AAE 
NAI 


20s 

25s 

30s 

0.10 

0.13 

0.16 

0.05 

0.07 

0.10 

0.07 

0.06 

0.08 
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5.4 


INVERSION  RESULTS 


The  dispersion  data  for  each  of  the  five  African  paths 
were  inverted  with  the  iterative  linear  inversion  technique 
employed  in  an  earlier  study  of  NTS  surface  waves  (Bache, 
et  al. , 1978).  With  this  technique,  phase  and  group  velocities 
are  simultaneously  fit  with  a smooth  plane-layered  model  of 
crustal  and  upper  mantle  shear  velocity  (8) , compressional  velo 
city  (a)  and  density  (p) . The  inversion  model  obtained  for  a 
given  path  represents  an  estimate  for  the  average  structure 
along  that  path. 

It  is  desirable  to  build  additional  constraints  into  the 
inversion  of  dispersion  data  in  order  to  supplement  their  re- 
solving power.  In  addition  to  the  smoothness  requirement,  our 
inversion  algorithm  allows  for  two  types  of  constraints.  First 
linear  a- 6 and  p-8  relationships  may  be  imposed  to  compensate 
for  the  fact  that  phase  and  group  velocities  are  only  weakly 
dependent  on  a and  p.  Second,  the  depth  to  discontinuities 
(such  as  the  crust-mantle  boundary)  may  be  incorporated  into 
the  model.  In  general,  dispersion  data  cannot  distinguish  a 
sharp  discontinuity  from  a gradual  variation  in  8 with  depth. 
These  were  the  only  constraints  used  to  invert  the  African 
data.  No  "guess  models"  were  used. 

The  model  constraints  employed  are  well-suited  for 
incorporating  results  from  previous  geophysical  surveys,  such 
as  Pn  velocities  and  crustal  thicknesses  determined  from 
refraction  studies.  While  some  information  of  this  type  is 
available  for  the  African  paths  studied  here,  it  was  not  suf- 
ficient for  confidently  fixing  the  constraints.  Instead,  we 
employed  "default"  constraints  whose  primary  purpose  was  to 
ensure  self-consistent  and  reasonable  inversion  models.  For 
each  of  the  five  paths  we  constrained  a and  p in  both  the 
crust  and  mantle  by 


a = 1.732  B (Poisson's  ratio  = 0.25) 

p = 0.3  a + 0.9  (a  form  of  Birch's  law).  (5.2) 

Using  Equation  (5.2),  the  data  for  each  path  were  inverted  for 
B with  and  without  assuming  a crustal  thickness.  For  some 
paths  two  or  more  crustal  thicknesses  were  tried  (30,  35  or 
40  km) . In  addition,  a 2 km  sedimentary  layer  was  incorporated 
into  the  inversion  models  for  the  HOG-HLW  and  HOG-X  paths  be- 
cause their  dispersion  data  seemed  to  require  low  shear  velo- 
cities near  the  surface. 

Figures  21  through  25  show  two  or  more  inversion  models 
obtained  for  each  of  the  five  African  paths.  The  different 
models  shown  for  each  path,  corresponding  to  different  assump- 
tions about  the  crustal  thickness,  predict  only  slightly  dif- 
ferent dispersion.  The  model  dispersion  is  compared  to  the 
observations  in  each  figure. 

With  the  exception  of  the  HOG-AAE  models,  the  inversion 
models  closely  agree  with  all  the  observed  velocity  points 
except  those  that  cannot  be  fit  by  a smooth  curve.  For  the 
HOG-AAE  path  the  inversion  models  do  not  simultaneously  fit 
the  Rayleigh  and  Love  wave  data.  Figure  23a  shows  three  inver- 
sion models  determined  from  only  the  AAE  Rayleigh  phase  and 
group  velocity  data.  They  match  these  data  well,  but  system- 
atically misfit  the  Love  wave  group  velocities.  An  inversion 
of  both  data  sets  together  (Figure  23b)  produces  a model  that 
fits  the  Love  wave  data  slightly  better,  but  the  Rayleigh  wave 
data  much  worse.  We  do  not  know  the  reason  for  this  incompati- 
bility between  the  Love  and  Rayleigh  data  from  AAE,  but  an 
apparent  anisotropy  due  to  lateral  variations  along  the  HOG- 
AAE  path  is  a plausible  possibility.  In  contrast,  a simulta- 
neous inversion  of  Rayleigh  and  Love  wave  data  from  the  HOG-HLW 
path  was  quite  successful. 
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Two  inversion  models  for  the  HOG-SHI  path  are  shown  (left)  together  with 
a comparison  of  the  observed  dispersion  with  the  predicted  dispersion  for 
the  model  for  which  crustal  thickness  was  not  assumed. 


Figure  23a.  Three  inversion  models  inverted  from  the  HOG-AAE  Rayleigh  wave  dispersion 
data  are  shown  at  left.  At  right,  the  predicted  dispersion  (Rayleigh  and 
Love)  for  the  model  with  a 35  km  crust  is  compared  to  the  observed 


Figure  23b.  A model  inverted  from  the  HOG-AAE  Rayleigh  and  Love  dispersion  data  with 
no  crustal  thickness  assumed  is  shown  (left)  together  with  a comparison 
of  its  predicted  dispersion  to  the  observed  dispersion. 


Figure  25.  Two  inversion  models  for  the  HOG-X  path  are  shown  (left)  together  with  a 
comparison  of  the  observed  dispersion  to  the  predicted  dispersion  for  the 
model  with  the  35  km  crust. 


We  found  that  models  with  either  no  sharp  crust-mantle 
transition  or  a crustal  thickness  of  35  km  could  fit  the  data 
for  each  path.  For  the  HLW  and  AAE  paths,  crustal  thicknesses 
of  30  and  40  km  also  fit  the  data  and  we  would  expect  this  to 
be  the  case  for  the  other  paths  as  well.  However,  different 
crustal  thicknesses  imply  different  shear  velocities  at  the 
base  of  the  crust  and  top  of  the  mantle  (e.g.,  Figures  21  and 
23a) . Therefore,  independent  information  about  lower  crustal 
velocity,  Sn  velocity,  or  velocity  gradients  in  the  crust  and 
mantle  would  allow  one  to  narrow  the  range  of  plausible  crustal 
thicknesses. 

The  results  of  previous  studies  of  African  structure 
provide  some  basis  for  judging  the  inversion  models  in  Figures 
21  through  25.  Gumper  and  Pomeroy  (1970)  analyzed  surface 
waves  and  short  period  body  waves  from  numerous  African  paths, 
with  many  in  eastern  Africa.  Their  inferred  crustal  and  upper 
mantle  model  (AFRIC)  is  listed  in  Table  6.  Gumper  and  Pomeroy 
reported  Lg  velocities  between  3.48  and  3.72  km/s.  Lg  velo- 
city is  an  estimate  for  average  crustal  shear  velocity  and  our 
models  are  consistent  with  the  wide  range  they  observed.  Their 
body-wave  observations  indicated  Sn  velocities  ranging  from 
4.55  to  4.72  km/s.  For  comparison,  the  approximate  Sn  velo- 
cities from  our  inversion  models  are  given  in  Table  7.  Our  Sn 
velocities  seem  to  be  low  for  some  of  the  paths.  However, 
Gumper  and  Pomeroy  noted  that  Sn  did  not  propagate  efficiently 
along  paths  crossing  the  African  rift,  near  which  AAE  and  NAI 
are  located.  Furthermore,  Knopoff  and  Schlue  (1972)  inferred 
low  Sn  velocities  (4.25  to  4.55  km/s)  for  a surface-wave  path 
from  AAE  to  NAI.  In  addition,  refraction  studies  near  AAE  by 
Searle  and  Gouin  (1971)  found  Sn  *»  4.29  km/s.  The  latter  study 
also  reported  a very  thick  crust,  up  to  48  km,  under  AAE.  A 
crust  this  thick  will  violate  our  HOG-AAE  data.  While  these 
studies  suggest  anomalous  structures  near  the  African  rift, 
we  must  keep  in  mind  that  only  a small  portion  of  the  HOG-AAE 
and  HOG-NAI  paths  traverses  the  rift  zone. 
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TABLE  6 

AFRIC  MODEL  OF  GUMPER  AND  POMEROY  (1970) 


Layer 

Thickness 

Depth  to 
Layer  Bottom 

a 

6 

P 

7.0 

7.0 

5.90 

3.35 

2.70 

10.5 

17.5 

6.15 

3.55 

2.80 

18.7 

36.2 

6.60 

3.72 

2.85 

80.0 

116.2 

8.05 

4.63 

3.30 

TABLE  7 

APPROXIMATE  Sn  VELOCITIES  FROM  INVERSION 
MODELS  WITH  35  KM  THICK  CRUST 


Path 

Sn  (kra/s) 

HOG-HLW 

4.6 

HOG-SHI 

4.6 

HOG-AAE 

in 

• 

HOG-NAI 

4.4 

HOG-X 

in 

• 
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While  the  inversion  models  found  for  each  path  are  rea- 
sonable, they  have  not  yet  been  integrated  to  give  a consistent 
picture  of  crustal  and  upper  mantle  structure  in  northern 
Africa.  For  example  even  though  the  HOG-AAE  and  HOG-NAI  paths 
are  fairly  close  (Figure  17) , the  inversion  models  for  these 
two  paths  are  more  different  than  we  would  expect.  However,  it 
may  be  that  very  similar  models  would  fit  both  sets  of  data 
just  as  well  considering  the  errors  in  the  data.  From  Figure 
17  we  also  note  that  the  HOG-HLW  path  is  essentially  the  first 
half  of  the  HOG-SHI  path.  This  fact  should  be  taken  into 
account. 

In  summary,  the  models  presented  in  this  section  repre- 
sent the  results  of  the  initial  inversion  of  the  dispersion 
data.  With  the  exception  of  the  Love-Rayleigh  inconsistency 
for  the  HOG-AAE  path  which  requires  more  study,  these  data  are 
fit  quite  well  by  smooth  and  reasonable  crust  and  upper  mantle 
models.  The  next  step  is  to  introduce  as  a kind  of  side  con- 
straint the  requirement  that  the  models  present  a coherent 
picture  of  the  structure  in  North  Africa.  This  will  lead  to 
our  final  models  for  the  paths  studied. 


5.5  SEISMOGRAM  SYNTHESIS 

At  this  point  in  our  study  of  surface  waves  from 
SAPHIRE,  synthetic  seismograms  generated  from  our  inversion 
models  serve  mainly  as  a check  on  our  data  analysis  and  inver- 
sion results,  and  as  a test  of  whether  the  observed  waveforms 
are  explainable  by  plane-layer  propagation  effects.  Before  we 
can  analyze  the  observed  amplitudes  and  infer  source  param- 
eters, we  must  resolve  certain  problems  with  propagation 
effects  and  obtain  a more  consistent  set  of  path  models. 

Figure  26  compares  the  observed  SAPHIRE  seismograms 
from  Figure  18  to  synthetic  surface  waves  generated  from  one 
of  the  inversion  models  for  each  path.  In  each  comparison. 


w 


HLW  Transverse  (A  = 2688  km) 


SHI  Vertical  (A  = 4729  km) 


Figure  26. 


seconds 


The  synthetic  seismograms  computed  from  the  inversion 
models  are  compared  to  the  observed  seismograms  from 
each  of  the  five  stations.  For  AAE  and  NAI,  the  syn- 
thetics obtained  from  two  y(<*))  models  (QAF  and  QAAE) 
are  shown. 
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the  observed  and  synthetic  seismograms  are  correctly  aligned 
in  time.  The  AAE  transverse  component  was  not  synthesized 
since  we  did  not  fit  the  AAE  Love  wave  dispersion  data. 

We  match  the  waveforms  of  the  HLW , SHI  and  X seismograms 
very  well.  (The  latter  part  of  the  observed  HLW  transverse 


r 


c 

TABLE 

8 

GAMMA  MODELS 

USED  IN  SURFACE-WAVE 

SYNTHESIS 

Yr  (10* 

/km) 

yl 

1 

Period 

QAF 

QAAE 

QAF 

40 

0.63 

4.00 

1.65 

] 1 _ 

1 1 c 

30 

0.591 

4.00 

1.30 

1 

25 

0.597 

4.00 

0.98 

20 

0.70 

3.15 

0.85 

1 C' 

18 

0.84 

2.59 

1.03 

16 

1.10 

1.96 

1.40 

14 

1.40 

1.25 

1.90 

c 

12 

1.80 

0.83 

2.90 

1 

10 

3.00 

2.20 

6.00 

8 

8.00 

4.44 

15.40 

o 

6 

12.00 

6.50 

22.20 

4 

20.00 

9.00 

36.00 

r 


We  are  not  now  certain  whether  reasonable  plane-layered 
velocity  and  attenuation  models  can  account  for  the  waveforms 
of  the  AAE  and  NAI  seismograms.  The  proximity  of  a major  tec- 
tonic feature,  the  African  rift,  may  play  an  important  role. 

In  summary,  the  synthetic  seismograms  computed  with  a 
simple  source  show  good  agreement  with  the  observations  at 
most  of  the  stations,  though  further  work  is  required.  As  we 
mentioned  at  the  end  of  Section  5.4,  we  first  need  to  complete 
our  development  of  consistent  structural  models  for  this 
region.  We  also  need  to  develop  reasonable  and  consistent  Q 
or  y models  to  account  for  the  attenuation.  Having  done  this 
we  will  use  the  comparison  between  observed  and  synthetic 
seismograms  to  infer  the  source  characteristics. 
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